Review: FT (Continuous) vs. DTFT (Discrete)

Why the Z-Transform?

• The Fourier Transform assumes signals are stable and well-behaved
• But some biosignals or systems may not be absolutely summable
• The Z-Transform generalizes the Fourier Transform
• Useful for analyzing discrete-time systems, especially when stability and causality matter

Definition

Let \(x[n]\) be a discrete-time signal.

The Z-Transform is defined as:

\[X(z) = \sum_{n=-\infty}^{\infty} x[n] z^{-n}\]

Where: - \(z \in \mathbb{C}\) is a complex variable - \(z = re^{j\omega}\)

Region of Convergence (ROC)

• The Z-Transform converges only for certain values of \(z\)
• The set of \(z\) for which the series converges is the ROC
• ROC is critical for system stability and causality

Z-Plane Representation

• Poles: values of \(z\) where \(X(z) \to \infty\)
• Zeros: values where \(X(z) = 0\)
• Visualization of poles and zeros helps in understanding system behavior

(-1.5, 1.5)

(-1.5, 1.5)

Relationship with Fourier Transform

If the ROC includes the unit circle, \(|z| = 1\), then:

\[X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}\]

So the Fourier Transform is a special case of the Z-Transform.

Example Signal

Let:

\[x[n] = a^n u[n]\]

Where: - \(a \in \mathbb{R}\) - \(u[n]\) is the unit step function (0 for \(n<0\), 1 for \(n\geq 0\))

Step 1: Z-Transform

Apply the definition:

\[X(z) = \sum_{n=0}^{\infty} a^n z^{-n} = \sum_{n=0}^{\infty} (az^{-1})^n\]

This is a geometric series:

\[X(z) = \frac{1}{1 - az^{-1}} = \frac{z}{z - a}\]

What is a Geometric Series?

A geometric series is a sum of terms where each term is multiplied by the same constant:

\[S = \sum_{n=0}^{\infty} r^n = 1 + r + r^2 + r^3 + \cdots\]

The value of \(r\) determines whether this sum converges (has a finite limit) or diverges.

Goal

Understand why this series:

\[\sum_{n=0}^{\infty} r^n\]

converges if and only if:

\[\boxed{|r| < 1}\]

Partial Sums

Let’s consider the sum up to the \(N\)-th term:

\[S_N = \sum_{n=0}^{N} r^n = 1 + r + r^2 + \cdots + r^N\]

This has a closed-form expression:

\[S_N = \frac{1 - r^{N+1}}{1 - r} \quad \text{for } r \neq 1\]

Taking the Limit

To find the sum of the infinite series, take the limit as \(N \to \infty\):

\[S = \lim_{N \to \infty} S_N = \lim_{N \to \infty} \frac{1 - r^{N+1}}{1 - r}\]

Case 1: \(|r| < 1\)

If \(|r| < 1\), then:

\[r^{N+1} \to 0 \quad \text{as } N \to \infty\]

So the sum becomes:

\[S = \frac{1}{1 - r}\]

✅ The geometric series converges.

Case 2: \(|r| \geq 1\)

• If \(r = 1\), then \(S_N = N + 1 \to \infty\)
• If \(r = -1\), the sum oscillates: \(1 - 1 + 1 - 1 + \cdots\)
• If \(|r| > 1\), then \(r^{N+1} \to \infty\)

❌ In all cases: the series diverges

Intuition

• When \(|r| < 1\), each term \(r^n\) gets smaller and smaller
• Their total sum settles to a finite number
• When \(|r| \geq 1\), the terms don’t vanish — the sum keeps growing or oscillating

Example: \(r = 0.5\)

\[S = 1 + 0.5 + 0.25 + 0.125 + \cdots = \frac{1}{1 - 0.5} = 2\]

Every term adds less. The sum “flattens out.”

Why This Matters

The Z-Transform often gives us geometric series like:

\[\sum_{n=0}^{\infty} (az^{-1})^n\]

This converges only if:

\[|az^{-1}| < 1 \Rightarrow |z| > |a|\]

So, understanding convergence of geometric series = understanding ROC in Z-transforms.

Summary

\[\sum_{n=0}^{\infty} r^n = \frac{1}{1 - r} \quad \text{if } |r| < 1\]

Step 2: Region of Convergence (ROC)

For convergence of the geometric series:

\[|az^{-1}| < 1 \Rightarrow |z| > |a|\]

Therefore, the ROC is:

\[\boxed{|z| > |a|}\]

Interpretation

• Causal signal (defined for \(n \geq 0\))
• ROC is outside the outermost pole
• Stable system only if ROC includes the unit circle \(|z| = 1\)

Example: \(a = 0.5\)

\[x[n] = (0.5)^n u[n]\]

Z-Transform:

\[X(z) = \sum_{n=0}^{\infty} (0.5)^n z^{-n} = \frac{z}{z - 0.5}\]

Region of Convergence:

\[\boxed{|z| > 0.5}\]

Summary

Properties of the Z-Transform

• Linearity: \(a x[n] + b y[n] \to aX(z) + bY(z)\)
• Time shifting: \(x[n - k] \to z^{-k} X(z)\)
• Scaling in the z-domain: \(a^n x[n] \to X(z/a)\)
• Convolution: \(x[n] * h[n] \to X(z)H(z)\)

Example

Let \(x[n] = a^n u[n]\), where \(|a| < 1\)

\[X(z) = \sum_{n=0}^{\infty} a^n z^{-n} = \frac{1}{1 - az^{-1}}, \quad \text{ROC: } |z| > |a|\]

Difference Equations in DSP

A difference equation relates input and output values at different time steps.

\[y[n] - a_1 y[n-1] - a_2 y[n-2] = b_0 x[n] + b_1 x[n-1]\]

Common in: - Digital filters (FIR, IIR) - Signal models in ECG, EEG analysis - Implementation in real-time biosignal systems

Z-Transform of Time-Shifted Terms

The Z-Transform turns time shifts into powers of \(z^{-1}\):

Step 1: Apply Z-Transform

Given:

\[y[n] - a_1 y[n-1] - a_2 y[n-2] = b_0 x[n] + b_1 x[n-1]\]

Apply \(\mathcal{Z} \{ \cdot \}\):

\[Y(z) - a_1 z^{-1} Y(z) - a_2 z^{-2} Y(z) = b_0 X(z) + b_1 z^{-1} X(z)\]

Step 2: Factor and Solve for \(H(z)\)

Group:

\[Y(z)(1 - a_1 z^{-1} - a_2 z^{-2}) = X(z)(b_0 + b_1 z^{-1})\]

Divide both sides:

\[H(z) = \frac{Y(z)}{X(z)} = \frac{b_0 + b_1 z^{-1}}{1 - a_1 z^{-1} - a_2 z^{-2}}\]

Example

Given:

\[y[n] - 0.9 y[n-1] = x[n] - 0.5 x[n-1]\]

Z-Transform:

\[Y(z)(1 - 0.9 z^{-1}) = X(z)(1 - 0.5 z^{-1})\]

Transfer Function:

\[H(z) = \frac{1 - 0.5 z^{-1}}{1 - 0.9 z^{-1}}\]

Poles and Zeros

Let’s analyze \(H(z)\):

• Zeros: Roots of the numerator \(\Rightarrow z = 0.5\)
• Poles: Roots of the denominator \(\Rightarrow z = 0.9\)

Practice

Convert this equation:

\[y[n] = 0.6 y[n-1] + x[n] + x[n-1]\]

Find: - \(H(z)\) - Poles and zeros - Plot them in the Z-plane

Application in Biosignal Processing

• Analysis of digital filters for ECG, EEG, etc.
• Design of stable and causal filtering systems
• Useful in difference equation modeling of biosignals

Summary

• Z-Transform is a powerful tool for analyzing discrete systems
• Provides insight into stability, causality, and system behavior
• A generalization of the Fourier Transform
• Crucial in digital signal processing of biosignals
• Z-Transform converts difference equations into algebraic expressions
• Transfer function \(H(z)\) tells us how the system responds to inputs
• Key for digital filter design in biosignal processing

Next Steps

• Practice Z-Transform computations
• Pole-zero plotting exercises
• Application to real biosignal filtering problems

Filter Design

Filter Design

def design_filter(zeros=None, poles=None, gain=1.0):
    """
    Diseña un filtro digital a partir de ceros y/o polos y una ganancia.

    Parámetros:
    - zeros: lista de ceros (raíces del numerador), o None para no incluir
    - poles: lista de polos (raíces del denominador), o None para no incluir
    - gain: ganancia escalar del filtro

    Devuelve:
    - b: coeficientes del numerador
    - a: coeficientes del denominador
    """
    # Si no se pasan ceros, asumimos un FIR trivial (b = [gain])
    if zeros:
        b = gain * np.poly(zeros)
    else:
        b = np.array([gain], dtype=float)

    # Si no se pasan polos, asumimos sistema FIR (a = [1])
    if poles:
        a = np.poly(poles)
    else:
        a = np.array([1.0], dtype=float)

    return b, a

Filter Design

def gaussian(x, mu, sigma, A):
    """
    Genera una función gaussiana.

    Parámetros:
    - x: array de tiempos
    - mu: posición central de la gaussiana
    - sigma: desviación estándar
    - A: amplitud
    """
    return A * np.exp(-((x - mu) ** 2) / (2 * sigma**2))

Filter Design

def simulate_ecg(duration=10.0, fs=500, heart_rate=60):
    """
    Simula un ECG sintético basado en la superposición de ondas gaussianas.

    Parámetros:
    - duration: duración de la señal en segundos
    - fs: frecuencia de muestreo en Hz
    - heart_rate: frecuencia cardiaca en latidos por minuto

    Devuelve:
    - t: vector de tiempos
    - ecg: señal simulada de ECG en mV
    """
    dt = 1 / fs
    t = np.arange(0, duration, dt)
    rr = 60 / heart_rate  # intervalo RR en segundos

    # Inicializar señal
    ecg = np.zeros_like(t)

    # Parámetros de las ondas (posiciones relativas en segundos)
    # P wave
    p_amp, p_dur, p_delay = 0.25, 0.09, 0.16
    # Q wave
    q_amp, q_dur, q_delay = -0.05, 0.066, 0.166
    # R wave
    r_amp, r_dur, r_delay = 1.6, 0.1, 0.166
    # S wave
    s_amp, s_dur, s_delay = -0.25, 0.066, 0.19
    # T wave
    t_amp, t_dur, t_delay = 0.35, 0.142, 0.36

    # Generar cada latido
    for beat_start in np.arange(0, duration, rr):
        mask = (t >= beat_start) & (t < beat_start + rr)
        tb = t[mask] - beat_start
        ecg[mask] += gaussian(tb, p_delay, p_dur / 2, p_amp)
        ecg[mask] += gaussian(tb, q_delay, q_dur / 2, q_amp)
        ecg[mask] += gaussian(tb, r_delay, r_dur / 2, r_amp)
        ecg[mask] += gaussian(tb, s_delay, s_dur / 2, s_amp)
        ecg[mask] += gaussian(tb, t_delay, t_dur / 2, t_amp)

    return t, ecg+0.1*np.cos(2*np.pi*60*t)

Filter Design

# Parámetros de simulación
DURATION = 10    # segundos
FS = 500         # Hz
HR = 70          # latidos por minuto

# Generar señal
t, ecg_signal = simulate_ecg(duration=DURATION, fs=FS, heart_rate=HR)

# Graficar resultado
plt.figure(figsize=(12, 4))
plt.plot(t, ecg_signal, linewidth=1)
plt.title(f'Señal de ECG sintética ({HR} bpm)')
plt.xlabel('Tiempo (s)')
plt.ylabel('Amplitud (mV)')
plt.grid(True)
plt.tight_layout()
plt.show()

Filter Design

n = len(ecg_signal)
yf = np.fft.fft(ecg_signal)
xf = np.fft.fftfreq(n, 1/fs)[:n//2]
plt.plot(xf, 2.0/n * np.abs(yf[0:n//2]))
plt.grid()
plt.show()

Filter Design

fc = 60
b,a =design_filter(zeros=[np.exp(1j*2*np.pi*fc/fs), np.exp(-1j*2*np.pi*fc/fs)])
print(a)
print(b)

w, h = sig.freqz(a, b, worN=8000)
plt.plot(w/np.pi*fs/2, 20*np.log10(abs(h)))
plt.grid()
plt.xlabel('Frequency (Hz)')

Filter Design

Filter Design

You obtain \(h_d[n]\) as the inverse discrete–time Fourier transform of your ideal frequency response \(H_d(e^{j\omega})\). For a low-pass filter with cutoff \(\omega_c\),

\[ H_d(e^{j\omega}) = \begin{cases} 1, & |\omega|\le\omega_c,\\ 0, & \text{otherwise}. \end{cases} \]

By definition of the inverse DTFT,

\[ h_d[n] \;=\; \frac{1}{2\pi}\int_{-\pi}^{\pi} H_d(e^{j\omega})\,e^{j\omega n}\,d\omega \;=\;\frac{1}{2\pi}\int_{-\omega_c}^{\omega_c} e^{j\omega n}\,d\omega. \]

Relating to sampling frequency

If your cutoff is specified in Hz, \(f_c\), and sampling rate is \(f_s\), then

\[ \omega_c = 2\pi\,\frac{f_c}{f_s}, \]

so you can write

\[ h_d[n] =\frac{\sin\!\bigl(2\pi \frac{f_c}{f_s}\,n\bigr)}{\pi\,n} =\;2\;\frac{f_c}{f_s}\;\frac{\sin\!\bigl(2\pi \frac{f_c}{f_s}\,n\bigr)}{2\pi \frac{f_c}{f_s}\,n} =2\frac{f_c}{f_s}\,\mathrm{sinc}\!\Bigl(2\frac{f_c}{f_s}\,n\Bigr), \]

where we define the normalized sinc as \(\mathrm{sinc}(x)=\frac{\sin(\pi x)}{\pi x}\).

Key takeaway

• \(h_d[n]\) is exactly the inverse‐DTFT of the ideal (“brick‐wall”) frequency specification.
• It produces a sinc-shaped impulse response of infinite length.
• Truncation (with a window) makes it realizable as a finite-length FIR.

Filter Design

The typical characteristics of classic windows are:

• Rectangular: narrowest main lobe (width ≈ 4π/M rad) but highest sidelobes (first sidelobe ≈ –13 dB, stop-band attenuation ≈ 21 dB). It gives the steepest transition for a given filter order, at the expense of poorer stop-band rejection.

• Bartlett (triangular): somewhat wider main lobe (≈ 8π/M), sidelobes ≈ –25 dB.

• Hann: main lobe ≈ 8π/M, sidelobes ≈ –31 dB (better rejection than rectangular but smoother transitions).

• Hamming: main lobe ≈ 8π/M, sidelobes ≈ –41 dB (minimum stop-band attenuation ≈ 53 dB). Very popular for its good compromise between transition width and stop-band attenuation.

• Blackman: wider main lobe (≈ 12π/M) but very low sidelobes (≈ –57 dB, attenuation ≈ 74 dB).

• Kaiser: allows selection of a parameter β to control sidelobe attenuation, offering flexibility. Approximately, to achieve A dB of attenuation,

  \[ \beta \approx 0.1102\,(A - 8.7)\quad(\text{for }A>50), \]

  and the normalized transition width Δω relates to the order M and β by

  \[ M \approx \frac{A - 8}{2.285\,\Delta\omega} \]

These formulas stem from Kaiser’s approximations and help in sizing the filter.

Windows Forms

(-100.0, 5.0)

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal.windows import (
    boxcar,       # Rectangular
    bartlett,    # Triangular
    hann,        # Hann
    hamming,     # Hamming
    blackman,    # Blackman
    kaiser       # Kaiser
)

# Parameters
M = 64               # window length
beta = 8.6           # Kaiser parameter for moderate sidelobe attenuation
nfft = 512           # FFT length for frequency response
fs = 1.0             # normalized sampling rate

# Generate windows
windows = {
    'Rectangular': boxcar(M),
    'Bartlett':    bartlett(M),
    'Hann':        hann(M),
    'Hamming':     hamming(M),
    'Blackman':    blackman(M),
    f'Kaiser (β={beta})': kaiser(M, beta)
}

# Time-domain plot
plt.figure(figsize=(8, 4))
for name, w in windows.items():
    plt.plot(np.arange(M), w, label=name)
plt.title('Window Functions — Time Domain')
plt.xlabel('Sample index n')
plt.ylabel('Amplitude')
plt.legend(loc='upper right')
plt.grid(True)
plt.tight_layout()

# Frequency-domain plot
plt.figure(figsize=(8, 4))
for name, w in windows.items():
    # Compute normalized frequency response
    W = np.fft.fft(w, nfft)
    W_mag = 20 * np.log10(np.abs(W) / np.max(np.abs(W)))
    freqs = np.fft.fftfreq(nfft, d=1/fs)
    # Only plot 0 ≤ f ≤ 0.5 (normalized Nyquist)
    idx = np.logical_and(freqs >= 0, freqs <= 0.5)
    plt.plot(freqs[idx], W_mag[idx], label=name)
plt.title('Window Functions — Frequency Response')
plt.xlabel('Normalized Frequency (cycles/sample)')
plt.ylabel('Magnitude (dB)')
plt.ylim(-100, 5)
plt.legend(loc='upper right')
plt.grid(True)
plt.tight_layout()

plt.show()

Filter Design

N = 101
fc = 30

# 2. Compute the ideal impulse response h_d[n] of a low-pass filter
n = np.arange(N)
M = (N - 1) / 2
# Using the normalized sinc: sinc(x) = sin(pi*x)/(pi*x)
h_d = (2 * fc / FS) * np.sinc(2 * fc * (n - M) / FS)

# 3. Choose a window w[n] (here: Hamming) and truncate
w = np.hamming(N)
h = h_d * w

# 4. Normalize to ensure unity gain at DC
h /= np.sum(h)

# 6. Filter it (only allowed library call)
y = sig.lfilter(h, 1.0, ecg_signal)

# 7. Plot input vs. output
plt.figure(figsize=(8, 4))
plt.plot(t, ecg_signal, label='Original signal')
plt.plot(t, y, label='Filtered signal', linewidth=2)
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Low-pass FIR (window method) – Cutoff 100 Hz')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

Step 0: The Analog Domain - Laplace Transform

• Context: Continuous-time signals and Linear Time-Invariant (LTI) systems.
• Tool: The Laplace Transform converts differential equations (time domain) to algebraic equations (s-domain).
• Transfer Function \(H(s)\): Represents the system in the s-domain. \[H(s) = \frac{Y(s)}{X(s)}\] where \(s = \sigma + j\Omega\) is the complex frequency (\(\Omega\) = analog angular frequency).

Step 0: The Analog Domain - s-Plane Analysis

• s-Plane: Complex plane where \(H(s)\) is analyzed.
• Poles & Zeros: Roots of the denominator & numerator of \(H(s)\), respectively. Their locations determine filter behavior.
• Frequency Response: Determined by \(H(j\Omega)\) (evaluating \(H(s)\) on the imaginary axis).
• Stability: For a causal system to be stable, all poles must be in the Left-Half Plane (LHP), i.e., \(\text{Re}\{s\} < 0\).

Step 0: The Analog Domain - Filter Prototypes

• Standardized, well-understood analog filter approximations:
  • Butterworth: Maximally flat passband.
  • Chebyshev Type I: Equiripple passband, monotonic stopband.
  • Chebyshev Type II: Monotonic passband, equiripple stopband.
  • Elliptic (Cauer): Equiripple in both passband and stopband (sharpest transition for a given order).
• Design usually starts with a normalized low-pass prototype.

The Bridge: Analog-to-Digital Mapping

• Core Task: Convert the analog filter \(H(s)\) (s-plane) into a digital filter \(H(z)\) (z-plane).
• Need a mapping that relates the continuous frequency \(\Omega\) to the discrete frequency \(\omega\).
• Crucially, the stable region (LHP in s-plane) should map to the stable region (inside the unit circle in z-plane).

Step 1: Digital Filter Specifications

• Define the desired characteristics of the digital filter:
  • Critical Frequencies:
    • Passband edge: \(\omega_p\)
    • Stopband edge: \(\omega_s\)
    • (Normalized digital angular frequencies, e.g., \(0 \le \omega \le \pi\))
  • Tolerances:
    • Max. passband ripple/attenuation: \(A_p\) (dB) or \(\delta_p\)
    • Min. stopband attenuation: \(A_s\) (dB) or \(\delta_s\)

Step 2: Method 1 - Impulse Invariance

• Concept: Match the impulse response: \(h_{digital}[n] \approx T \cdot h_{analog}(nT)\).
• Mapping: Maps s-plane pole \(s_k\) to z-plane pole \(z_k = e^{s_k T}\).
  • Uses partial fraction expansion of \(H(s)\).
• Pros: Preserves impulse response shape.
• Cons: Suffers from frequency aliasing if \(H(j\Omega)\) is not bandlimited below Nyquist frequency (\(F_s/2\)). Best for low-pass/band-pass filters.

Step 2: Method 2 - Bilinear Transformation

• Concept: An algebraic substitution mapping the \(j\Omega\) axis onto the unit circle.
• Mapping Formula: \[s = \frac{2}{T} \frac{1 - z^{-1}}{1 + z^{-1}}\] (T = sampling period, often set to 2 for simplicity during derivation).
• Pros:
  • No aliasing: Maps entire \(j\Omega\) axis to the unit circle (\(-\pi \le \omega \le \pi\)).
  • Preserves stability: Maps LHP (stable s-region) to inside the unit circle (stable z-region).
• Cons: Introduces non-linear frequency warping.

Understanding Frequency Warping (Bilinear)

• The relationship between analog (\(\Omega\)) and digital (\(\omega\)) frequencies is non-linear: \[\Omega = \frac{2}{T} \tan\left(\frac{\omega T}{2}\right)\] or equivalently: \[\omega = 2 \arctan\left(\frac{\Omega T}{2}\right)\]
• This compresses the infinite analog frequency axis onto the finite digital frequency range \([-\pi, \pi]\).
• The mapping is non-uniform, most distorted near \(\omega = \pm \pi\).

Step 3: Frequency Pre-warping (Bilinear Only)

• Problem: Frequency warping distorts the locations of \(\omega_p\) and \(\omega_s\).
• Solution: Pre-warp the digital specifications (\(\omega_p, \omega_s\)) into corresponding analog frequencies (\(\Omega_p, \Omega_s\)) before designing the analog filter.
• Formula: Use the warping relationship: \[\Omega_p = \frac{2}{T} \tan\left(\frac{\omega_p T}{2}\right)\] \[\Omega_s = \frac{2}{T} \tan\left(\frac{\omega_s T}{2}\right)\]
• Use these \(\Omega_p, \Omega_s\) values for the analog design step.

1. Define Filter Specifications

This is the most crucial first step.

• Type: Low-Pass (LPF), High-Pass (HPF), Band-Pass (BPF), Band-Stop (BSF).
• Critical Frequencies (Analog, \(\Omega\) in rad/s or \(f\) in Hz):
  • Cutoff (\(\Omega_c, f_c\))
  • Passband Edge(s) (\(\Omega_p, f_p\))
  • Stopband Edge(s) (\(\Omega_s, f_s\))
• Attenuation Levels (dB):
  • Max. Passband Attenuation/Ripple (\(A_{max}\) or gpass)
  • Min. Stopband Attenuation (\(A_{min}\) or gstop)
• Phase Response: Linear phase (e.g., Bessel) or not critical.
• Impedance: Source/Load matching.

1. Specifications: Python Example (Low-Pass)

Let’s define specs for an analog low-pass filter. Frequencies for analog filters in scipy are typically in rad/s.

Passband edge: 6283.19 rad/s (1000 Hz)

Stopband edge: 9424.78 rad/s (1500 Hz)

Max Passband Loss: 1.0 dB

Min Stopband Attenuation: 40.0 dB

2. Choose an Approximation Method (Filter Family)

Select based on trade-offs (roll-off, ripple, phase).

• Butterworth: Maximally flat passband, smooth roll-off.
• Chebyshev Type I: Steeper roll-off than Butterworth, passband ripple.
• Chebyshev Type II: Steeper roll-off, stopband ripple, flat passband.
• Elliptic (Cauer): Steepest roll-off, ripple in both bands.
• Bessel-Thomson: Best phase linearity (constant group delay), slowest roll-off.

For our Python examples, we’ll use Butterworth.

Butterworth Filter Fundamentals

The Butterworth filter is maximally flat in the passband.

The normalized prototype has a cutoff frequency \(\omega_c = 1\) rad/s.

Poles are placed on the unit circle in the \(s\)-plane at angles: \[s_k = e^{j\pi\frac{2k + n - 1}{2n}}, \quad k = 1,2,\ldots,n.\]

The Butterworth polynomial of order \(n\) is: \[B_n(s) = \prod_{k=1}^{n} (s - s_k).\]

Normalized Lowpass (LP) Transfer Function

The normalized \(n\)th-order lowpass prototype is: \[H_{LP}(s) = \frac{1}{B_n(s)} = \frac{1}{\prod_{k=1}^n (s - s_k)}.\] This yields a cutoff at \(|H_{LP}(j1)| = 1/\sqrt{2}\) (the –3 dB point).

Normalized Highpass (HP) Transfer Function

Apply the frequency transformation \(s \to 1/s\) and multiply by \(s^n\): \[H_{HP}(s) = H_{LP}(1/s)\,s^n = \frac{s^n}{B_n(1/s)} = \frac{s^n}{\prod_{k=1}^n (1/s - s_k)} = \frac{\prod_{k=1}^n(-s_k)\;s^n}{\prod_{k=1}^n(s - s_k^{-1})}.\] This prototype has a normalized cutoff at \(\omega=1\) rad/s in the highpass sense.

Normalized Bandpass (BP) Transfer Function

Use the bandpass transformation: \[s \to \frac{s^2 + \Omega_0^2}{B\,s},\] where \(\Omega_0 = \sqrt{\omega_1\omega_2}\) and \(B = \omega_2 - \omega_1\).

For a normalized prototype, we set \(\Omega_0 = 1\), \(B = 1\), giving: \[H_{BP}(s) = H_{LP}\left(\frac{s^2 + 1}{s}\right) = \frac{(s)^n}{\prod_{k=1}^n\left(\frac{s^2 + 1}{s} - s_k\right)} = \frac{s^n}{\prod_{k=1}^n\left(s^2 + 1 - s\,s_k\right)}.\]

Normalized Bandstop (BS) Transfer Function

Use the bandstop transformation: \[s \to \frac{s}{\frac{s^2 + 1}{s}} = \frac{s^2}{s^2 + 1},\] or equivalently \(s \to (B\,s)/(s^2 + \Omega_0^2)\) with \(B=1\), \(\Omega_0=1\). Then: \[H_{BS}(s) = H_{LP}\left(\frac{s}{s^2 + 1}\right) = \frac{1}{B_n\left(\frac{s}{s^2+1}\right)} = \frac{1}{\prod_{k=1}^n\left(\frac{s}{s^2+1} - s_k\right)} = \frac{(s^2+1)^n}{\prod_{k=1}^n\left(s - s_k(s^2+1)\right)} = \frac{(s^2+1)^n}{\prod_{k=1}^n\left(s - s_ks^2 - s_k\right)}.\]

3. Determine the Filter Order (N)

Estimating the filter order (\(N\)) is a crucial first step in analog filter design after defining specifications. The order determines the filter’s complexity and steepness of its frequency response roll-off.

These equations are for Butterworth filters. The calculated \(N\) is a real number and must be rounded up to the next integer.

Common Variables: * \(N\): Filter order (integer, result of formula is rounded up). * \(A_{min}\): Minimum stopband attenuation (dB, positive value, e.g., 40 dB). * \(A_{max}\): Maximum passband attenuation/ripple (dB, positive value, e.g., 1 dB). * \(\log_{10}\): Base-10 logarithm. * All \(\Omega\) values are angular frequencies in rad/s.

Butterworth Low-Pass (LP) Filter Order

\[N \ge \frac{\log_{10}\left(\frac{10^{A_{min}/10}-1}{10^{A_{max}/10}-1}\right)}{2 \log_{10}\left(\frac{\Omega_s}{\Omega_p}\right)}\]

Explanation of LP-Specific Variables: * \(\Omega_p\): Passband edge angular frequency. The frequency up to which signals pass with at most \(A_{max}\) attenuation. * \(\Omega_s\): Stopband edge angular frequency. The frequency from which signals are attenuated by at least \(A_{min}\). * For a low-pass filter, \(\Omega_s > \Omega_p\). The term \(\frac{\Omega_s}{\Omega_p}\) is the transition ratio.

Butterworth High-Pass (HP) Filter Order

\[N \ge \frac{\log_{10}\left(\frac{10^{A_{min}/10}-1}{10^{A_{max}/10}-1}\right)}{2 \log_{10}\left(\frac{\Omega_p}{\Omega_s}\right)}\]

Explanation of HP-Specific Variables: * \(\Omega_p\): Passband edge angular frequency. The frequency from which signals pass with at most \(A_{max}\) attenuation. * \(\Omega_s\): Stopband edge angular frequency. The frequency down to which signals are attenuated by at least \(A_{min}\). * For a high-pass filter, \(\Omega_p > \Omega_s\). The term \(\frac{\Omega_p}{\Omega_s}\) is the transition ratio (inverted compared to LP to keep it >1 for the logarithm).

Butterworth Band-Pass (BP) Filter Order

The order \(N_{BP}\) is the same as an equivalent Low-Pass Prototype. \[N_{BP} \ge \frac{\log_{10}\left(\frac{10^{A_{min}/10}-1}{10^{A_{max}/10}-1}\right)}{2 \log_{10}\left(\Omega_{s,LP_equiv}\right)}\]

Explanation of BP-Specific Variables: * \(\Omega_{p1}, \Omega_{p2}\): Passband edge angular frequencies (\(\Omega_{p1} < \Omega_{p2}\)). * \(\Omega_{s1}, \Omega_{s2}\): Stopband edge angular frequencies (\(\Omega_{s1} < \Omega_{p1}\) and \(\Omega_{s2} > \Omega_{p2}\)). * \(\Omega_0 = \sqrt{\Omega_{p1}\Omega_{p2}}\): Geometric center frequency of the passband. * \(B = \Omega_{p2} - \Omega_{p1}\): Bandwidth of the passband. * \(\Omega_{s,LP_equiv} = \min\left( \left| \frac{\Omega_{s1}^2 - \Omega_0^2}{\Omega_{s1} B} \right|, \left| \frac{\Omega_{s2}^2 - \Omega_0^2}{\Omega_{s2} B} \right| \right)\) * This \(\Omega_{s,LP_equiv}\) is the stopband edge of an equivalent low-pass prototype whose passband edge is normalized to 1 rad/s. It represents the more stringent of the two possible transitions from the passband edges to the stopband edges.

Butterworth Band-Stop (BS) Filter Order

The order \(N_{BS}\) is the same as an equivalent Low-Pass Prototype. \[N_{BS} \ge \frac{\log_{10}\left(\frac{10^{A_{min}/10}-1}{10^{A_{max}/10}-1}\right)}{2 \log_{10}\left(\Omega_{trans,LP_equiv}\right)}\]

Explanation of BS-Specific Variables: * \(\Omega_{s1}, \Omega_{s2}\): Stopband edge angular frequencies defining the rejected band (\(\Omega_{s1} < \Omega_{s2}\)). \(A_{min}\) is the attenuation in this band. * \(\Omega_{p1}, \Omega_{p2}\): Passband edge angular frequencies outside the stopband (\(\Omega_{p1} < \Omega_{s1}\) and \(\Omega_{p2} > \Omega_{s2}\)). \(A_{max}\) is the ripple in these passbands. * \(\Omega_0 = \sqrt{\Omega_{s1}\Omega_{s2}}\): Geometric center frequency of the stopband. * \(B = \Omega_{s2} - \Omega_{s1}\): Bandwidth of the stopband. * \(\Omega_{p,LP_i} = \left| \frac{\Omega_{pi} B}{\Omega_{pi}^2 - \Omega_0^2} \right|\) for \(i=1\) (using \(\Omega_{p1}\)) and \(i=2\) (using \(\Omega_{p2}\)). * These are the passband edges of an equivalent low-pass prototype whose stopband edge is normalized to 1 rad/s. * \(\Omega_{trans,LP_equiv} = \frac{1}{\min(\Omega_{p,LP1}, \Omega_{p,LP2})}\) * This is the effective transition ratio (stopband edge / passband edge) for the equivalent low-pass prototype.

3. Determine the Filter Order (N) & Denormalization

Calculated based on specs and chosen filter type. scipy.signal.buttord (and similar for other types) can find this for analog filters.

Required Butterworth Filter Order (N): 14

Butterworth Natural Frequency (Wn_buttord): 6593.83 rad/s

• Wn_buttord is the natural angular frequency (\(\Omega_n\)) for the filter. For Butterworth, this is the -3dB cutoff frequency \(\Omega_c\).

3. Determine the Filter Order (N) & Denormalization

Normalized Lowpass to Denormalized Lowpass

• Goal: Transform a normalized LP (\(\omega_c = 1\)) to a denormalized LP with cutoff \(\Omega_c\).
• Transformation: Replace \(s_{norm}\) with \(s / \Omega_c\). \[s_{norm} \to \frac{s}{\Omega_c}\]
• Denormalized Transfer Function: \[H_{LP}(s) = H_{norm}\left(\frac{s}{\Omega_c}\right)\]

Normalized Lowpass to Denormalized Highpass

• Goal: Transform a normalized LP (\(\omega_c = 1\)) to a denormalized HP with cutoff \(\Omega_c\).
• Transformation: Replace \(s_{norm}\) with \(\Omega_c / s\). \[s_{norm} \to \frac{\Omega_c}{s}\]
• Denormalized Transfer Function: \[H_{HP}(s) = H_{norm}\left(\frac{\Omega_c}{s}\right)\]

Normalized Lowpass to Denormalized Bandpass

• Goal: Transform a normalized LP (\(\omega_c = 1\)) to a denormalized BP with center frequency \(\Omega_0\) and bandwidth \(B\).
• Transformation: Replace \(s_{norm}\) with \(\frac{s^2 + \Omega_0^2}{B s}\). \[s_{norm} \to \frac{s^2 + \Omega_0^2}{B s}\]
• Parameters:
  • \(\Omega_0 = \sqrt{\Omega_1 \Omega_2}\) (geometric mean)
  • \(B = \Omega_2 - \Omega_1\) (bandwidth) (\(\Omega_1, \Omega_2\) are the desired band edges in rad/s)
• Denormalized Transfer Function: \[H_{BP}(s) = H_{norm}\left(\frac{s^2 + \Omega_0^2}{B s}\right)\]

Normalized Lowpass to Denormalized Bandstop

• Goal: Transform a normalized LP (\(\omega_c = 1\)) to a denormalized BS with center frequency \(\Omega_0\) and bandwidth \(B\).
• Transformation: Replace \(s_{norm}\) with \(\frac{B s}{s^2 + \Omega_0^2}\). \[s_{norm} \to \frac{B s}{s^2 + \Omega_0^2}\]
• Parameters:
  • \(\Omega_0 = \sqrt{\Omega_1 \Omega_2}\)
  • \(B = \Omega_2 - \Omega_1\)
• Denormalized Transfer Function: \[H_{BS}(s) = H_{norm}\left(\frac{B s}{s^2 + \Omega_0^2}\right)\]

4 & 5. Obtain Analog Filter Transfer Function \(H_a(s)\)

Conceptually: Find normalized LP prototype \(\rightarrow\) Denormalize \(\rightarrow\) Transform.

scipy.signal functions like butter, cheby1, etc., (with analog=True) perform these steps internally to give the final analog filter transfer function \(H_a(s)\).

\(H_a(s)\) can be represented as: 1. Numerator (\(b\)) and Denominator (\(a\)) coefficients: \(H_a(s) = \frac{b_0 s^M + b_1 s^{M-1} + \dots + b_M}{a_0 s^N + a_1 s^{N-1} + \dots + a_N}\) 2. Zeros (\(z\)), Poles (\(p\)), and Gain (\(k\)): \(H_a(s) = k \frac{\prod (s-z_i)}{\prod (s-p_j)}\)

4 & 5. Python: Design \(H_a(s)\) (Butterworth LPF)

Using signal.butter with analog=True. The Wn parameter here should be the cutoff frequency \(\Omega_c\). For buttord, the returned Wn_buttord is exactly this \(\Omega_c\).

Analog Filter Coefficients (Numerator b, Denominator a):

b_analog = [2.93718174e+53]

a_analog = [1.00000000e+00 5.88921844e+04 1.73414469e+09 3.37531810e+13
 4.84169654e+17 5.40639098e+21 4.84125527e+25 3.52958830e+29
 2.10491137e+33 1.02201934e+37 3.97946840e+40 1.20619642e+44
 2.69441500e+47 3.97843852e+50 2.93718174e+53]

Analog Filter Zeros, Poles, Gain (ZPK):

Zeros (z_analog) = []

Poles (p_analog) = [ -738.27500968+6552.37193897j -2177.8048373 +6223.80864963j
 -3508.13043664+5583.15760623j -4662.54372722+4662.54372722j
 -5583.15760623+3508.13043664j -6223.80864963+2177.8048373j
 -6552.37193897 +738.27500968j -6552.37193897 -738.27500968j
 -6223.80864963-2177.8048373j  -5583.15760623-3508.13043664j
 -4662.54372722-4662.54372722j -3508.13043664-5583.15760623j
 -2177.8048373 -6223.80864963j  -738.27500968-6552.37193897j]

Gain (k_analog) = 293718173729774468825490597908032846661500110135885824.00

Visualizing Analog Frequency Response (Python)

Use scipy.signal.freqs to compute the frequency response of an analog filter given its \(b,a\) coefficients.

(100.0, 7500)

(-60.0, 5.0)

(100.0, 7500)

6. Synthesize the Circuit (Realization)

Convert the mathematical \(H_a(s)\) into a physical electronic circuit. This step is not done by scipy.signal. It requires circuit design knowledge.

• Passive Filters:
  • Resistors (R), Inductors (L), Capacitors (C).
  • Commonly “ladder” structures.
  • Inductors can be bulky/non-ideal, especially at low frequencies.
• Active Filters:
  • Op-Amps, Resistors (R), Capacitors (C) (avoids inductors).
  • Can provide gain, easy to cascade.
  • Requires power. Limited by op-amp performance.
  • Common topologies: Sallen-Key, Multiple Feedback (MFB), State-Variable/Biquad.
  • Higher-order filters are often cascaded 2nd-order sections.

7. Component Selection and Simulation

• Component Values: Choose standard, available component values (R, C, L, Op-Amps) that are close to calculated ideals.
• Tolerances: Account for component variations.
• Non-Idealities: Consider real-world behavior of components.
• Circuit Simulation (Essential):
  • Use SPICE-based simulators (LTspice, QUCS, PSpice, etc.).
  • Verify performance against specifications.
  • Iterate on design and component values as needed.

This is also outside the scope of scipy.signal but critical for hardware implementation.

Tools for Analog Filter Design

• Mathematical Software:
  • Python with scipy.signal: Excellent for determining filter parameters (\(N, \Omega_c\)) and transfer functions (\(b,a\) or \(z,p,k\)), and plotting responses.
  • MATLAB (Signal Processing Toolbox), Octave.
• Filter Design Calculators/Software:
  • Texas Instruments FilterPro™, Analog Devices Filter Wizard®.
  • Online calculators and dedicated filter design suites.
• Circuit Simulators:
  • LTspice, QUCS, PSpice, NI Multisim™, Keysight ADS.
• Filter Design Handbooks:
  • Classic texts (e.g., by Van Valkenburg, Zverev) provide theory, tables, and formulas.

Step 5: Apply Digital Transformation

• Convert the designed analog filter \(H_a(s)\) to the digital domain \(H(z)\):
  • If Impulse Invariance: Use partial fractions and the \(s_k \rightarrow e^{s_k T}\) mapping.
  • If Bilinear Transformation: Substitute \(s = \frac{2}{T} \frac{1 - z^{-1}}{1 + z^{-1}}\) into \(H_a(s)\) and simplify algebraically.
• Result: The digital filter transfer function: \[H(z) = \frac{\sum_{k=0}^{M} b_k z^{-k}}{1 + \sum_{k=1}^{N} a_k z^{-k}}\]

Step 6: Realization (Implementation)

• The transfer function \(H(z)\) corresponds directly to a difference equation: \[y[n] = \sum_{k=0}^{M} b_k x[n-k] - \sum_{k=1}^{N} a_k y[n-k]\]
• \(x[n]\): Input sequence
• \(y[n]\): Output sequence
• \(b_k, a_k\): Filter coefficients derived from \(H(z)\).
• The feedback term (sum involving \(y[n-k]\)) makes the impulse response potentially infinite (IIR).

Step 1: Digital Filter Specifications

• Define the desired characteristics of the digital filter:
  • Sampling Frequency: \(F_s\) (Hz)
  • Critical Frequencies (Digital):
    • Passband edge: \(f_p\) (Hz) \(\implies \omega_p = 2\pi f_p / F_s\) (radians/sample)
    • Stopband edge: \(f_s\) (Hz) \(\implies \omega_s = 2\pi f_s / F_s\) (radians/sample)
    • Note: Scipy often uses frequencies normalized to Nyquist: \(W_p = f_p / (F_s/2)\), \(W_s = f_s / (F_s/2)\)
  • Tolerances:
    • Max. passband ripple/attenuation: \(A_p\) or gpass (dB)
    • Min. stopband attenuation: \(A_s\) or gstop (dB)

Step 1: Example Specifications (Python)

Sampling Frequency: 10000 Hz

Normalized Passband Edge: 0.200 (pi rad/sample)

Normalized Stopband Edge: 0.300 (pi rad/sample)

Max Passband Loss: 1 dB

Min Stopband Attenuation: 40 dB

Steps 2-4: Finding Order & Designing (Conceptual)

• Transformation: Bilinear Transform is implicitly used by many scipy.signal functions for IIR design.
• Frequency Pre-warping: Handled internally by Scipy when given digital frequency specs.
• Analog Design: Scipy calculates the required analog prototype order and parameters based on the (internally pre-warped) specs.
• Goal: Find filter order \(N\) and digital cutoff frequency \(W_n\).

Steps 2-4: Finding Order & Wn (Python - Butterworth Example)

• Use scipy.signal.buttord to find the minimum order and cutoff frequency for a digital Butterworth filter meeting the specs.

Minimum Filter Order (N): 12

Butterworth Natural Frequency (Wn): 0.211 (pi rad/sample)

• analog=False specifies we want parameters for a digital filter.
• Wn is the digital cutoff frequency (scalar for LP/HP, tuple for BP/BS) required by the butter function. It’s often close to wp.

Step 5: Designing the Digital Filter (Python - Butterworth)

• Use scipy.signal.butter (or cheby1, cheby2, ellip) to get the filter coefficients \(b\) (numerator) and \(a\) (denominator) of \(H(z)\).

Filter Coefficients:

Numerator (b): [0.0e+00 0.0e+00 1.0e-05 4.0e-05 1.0e-04 1.6e-04 1.9e-04 1.6e-04 1.0e-04
 4.0e-05 1.0e-05 0.0e+00 0.0e+00]

Denominator (a): [ 1.000000e+00 -6.930840e+00  2.271843e+01 -4.632962e+01  6.523097e+01
 -6.662310e+01  5.050575e+01 -2.858485e+01  1.197050e+01 -3.612930e+00
  7.452400e-01 -9.424000e-02  5.520000e-03]

• output='ba' gives numerator/denominator coefficients. Other options: 'sos' (second-order sections - numerically better), 'zpk' (zeros, poles, gain).

Visualizing Frequency Response (Python)

• Calculate the frequency response using scipy.signal.freqz.
• Plot magnitude (in dB) and phase.

(0.0, 5000.0)

(-60.0, 5.0)

(0.0, 5000.0)

Step 6: Realization & Filtering (Python Example)

• The coefficients b, a define the difference equation.
• Use scipy.signal.lfilter to apply the filter to a signal.

(-1.5, 1.5)